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Abstract. We describe, for the first time, a completely rigorous homotopy (path-following) algo- 
rithm (in the Turing machine model) to find approximate zeros of systems of polynomial equations. 
If the coordinates of the input systems and the initial zero are rational our algorithm involves only 
rational computations and if the homotopy is well posed an approximate zero with integer coordi- 
nates of the target system is obtained. The total bit complexity is linear in the length of the path 
in the condition metric, and polynomial in the logarithm of the maximum of the condition number 
along the path, and in the size of the input. 



1. Introduction 

The research on solving systems of polynomial equations has experienced a rapid development 
in the last two decades, both from a theoretical and from a practical perspective. Many of the 
recent advances are based on the study of the general idea of homotopy continuation methods: 
let / be the system whose solutions we want to find, and let g be another system whose solutions 
we already know. Then, join g and / with a homotopy, that is a curve ft in the vector space of 
polynomial systems, such that fo = g and /i = /, and try to follow the curves (homotopy paths) 
produced as a solution C,q of g is continued along the homotopy to a solution Q of ft- When t 
approaches 1, an approximation of a zero Ci of / is obtained. 

In order to describe such a method explicitly, we need two essential ingredients: 

(1) A construction of the starting system g and the homotopy path ft-, and, 

(2) once this path ft is chosen, a procedure to approximate Q (for a finite sequence of values 
of t starting with t = and ending with t = 1). 

The first of these two ingredients has been intensively studied from many perspectives, mainly 
using hnear homotopy paths, i.e. once ((y',Co) is chosen, C,q a solution of g, we just consider the 
path /t = (1 — t)g + tf. In |l2] a particularly simple choice of {gXo) was conjectured to be a 
good candidate for a initial pair, i.e. the complexity of homotopy methods with this starting pair 
could be polynomial on the average. This is still a challenging open conjecture that has been 
experimentally confirmed in [5]. In [6], [7], [8] it was proved that randomly chosen pairs ((?, Co) 
guarantee average polynomial complexity. In [13] a system whose zeros have coordinates equal 
to the roots of unity of appropriate degrees was proved to guarantee average quasi-polynomial 
complexity (polynomial for fixed degree.) 

In this paper we deal with the second of the two questions above, restricting ourselves to the 
case when the homotopy path that is followed is regular, i.e., Q is a regular zero of ft for every 
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Figure 1. Path-jumping scenario: we start following the path corresponding to a 
certain solution of /o but we may jump to another path in the middle. 



There exist several software packages which perform the path-following task of item ^ above 
(here is an incomplete hst: Bertini [1], HOM4PS2 |2lj, NAG4M2 |27|, and PHCpack |16]). In 
general, an initial step to is chosen, and a predictor step (a numerical integration step of the 
differential equation d{ft{Ct))/dt = 0) is made to approximate (to- A corrector step (several steps 
of Newton's method) is then used to get a better approximation of (to- This process is repeated 
by choosing ti,t2,... until t = 1 is reached. The software implementations mentioned above 
achieve spectacular practical results, with huge systems solved in a surprisingly short time. As 
a drawback, these fast methods make heuristic choices (notably the choice of ti), which may 
introduce uncertainty in the quality of the solutions they provide: how close to an actual zero is 
the given output? is the method actually following the path Q or maybe a path-jumping occurred 
in the middle? 

In Figure [T] we illustrate a path-jumping phenomenon that may occur when a heuristic predictor- 
corrector path-tracking procedure is used. 

In the problems with aim at computing all target solutions, this scenario can be detected simply 
by observing that two approximation sequences produce approximations to the same regular zero. 
The Kim-Smale a-test from [26], can be used to make this task rigorously, see [20] for an 
implementation of that test. Once detected, this can be remedied by rerunning the heuristic 
procedure with tighter tolerances and higher precision of the computation. 

However, an analysis of the end solutions would fail to detect the shortcomings of a heuristic 
method in the scenario with two approximation sequences "swapping" two paths as in Figure [2] 

In [12], a method which guarantees that path-jumping does not occur was shown for the first 
time, and its complexity (number of homotopy steps) was bounded above by a quantity depending 
on the maximum of the so-called condition number along the path {ft,(t), see (2.2) below. This 
result was recently improved in [39j, changing the maximum of the condition number along the path 
to the length of the path {ft,(t) in the "condition metric" (see Section [Z4 for details). However, 
the result in [39] does not fully describe an algorithm, for the explicit choice of the steps ti is not 
given. Describing a way to actually choose these steps is a nontrivial task that can be done in 
several fashions. There are three independent papers doing this job: [1], [T3], [17]. We briefly 
summarize in table [T] the properties of the algorithms in those papers and of the algorithm in this 
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Figure 2. Path-swapping scenario: two path-jumps occur producing correct (al- 
though permuted) target solutions. 

Table 1. The existing methods with certified output and complexity analysis for 
path-tracking. By arithmetic on M we mean that the BSS computation model [12] 
is assumed, that is exact arithmetical operations between real numbers are allowed. 
CL. means length in the condition metric. 



Method 


Complexity 


Arithmetic 


Assumptions on ft and comments 


[39J 


CL. 


R 


paths ft- Not constructive. 




CL. 


R 




US] 


>^C.L. 


R 


Linear paths. 


[T7] 


CL. 


R 


paths ft with special properties. 


This paper 


CL. 


Q 


Linear paths. 



paper too. The proof of the algorithm in [131 probably the shortest of the ones mentioned. As 
a small drawback, its complexity is not bounded above by the length of the path (/t,Ct) in the 
condition metric, but by the integral of the square of the condition number along the path ft, 
which is in general a (non-dramatically) greater quantity. 

All methods in Table [T] are originally designed for systems of homogeneous polynomials, and 
then an argument like that in [7j is used to produce affine approximate zeros of the original system, 
if this one was not homogeneous. 

Another difference from the heuristic methods described above is that there is no predictor step 
in these methods; only Newton's method is used (more exactly, projective Newton's method |38] 
described below): once to is chosen, one obtains Zt^ as the result of (projective) Newton's method 
with base system ft^ and base point zq = (q (or zq an approximate zero of Co)- The idea is repeated, 
again, until t = 1 is reached. 

Using the algorithm of [1], the main result of |39] reads as follows. 

Theorem 1. [39], [1] Assuming exact numerical computations, and assuming that the path ft is a 
great circle in the sphere § in the space of homogeneous systems, the steps tQ,ti, . . . can be chosen 
in such a way that the homotopy method outlined above produces an approximate zero of f = fi 
with associated exact zero Ci- The total number of steps is at most a small constant IXS'I'^ (d the 
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maximum of the degrees of the polynomials in f) times Cq = Colft, (t), that is the length of the path 
iftXt) in the so-called condition metric (if that length is infinity, the algorithm may never finish.) 

Here, the sphere § is the set of systems / such that ||/|| = 1 where || ■ || is the Bombieri-Weyl 
norm described below, and the condition metric is the usual product metric in § x P(C"'"'"^), point- 



wise multiplied by the condition number squared, see (2.3) below for a precise definition. The main 
result of ^39], [1] also applies to paths which are not great circles, although the constant Tld^^"^ 
may then vary. 

In recent works [5] , [27] we described a practical implementation of the "certified" method of jl] 
and used it in various experiments for estimating the average complexity of homotopy tracking. 

As we have already pointed out. Theorem [l] (and the other existing algorithms for path-tracking) 
needs exact numerical computations on real numbers. More precisely, it assumes computations 
under the BSS model of computation, [12], [TT]. While this assumption permits simplified analysis 
of the method, it is unrealistic in practice. Our main result here is to remove this assumption: if 
the input polynomials g and / and the initial approximate zero zq are given in (Gaussian) rational 
coordinates, then all the computations can be made over the rationals. We center our attention in 
linear paths, that is paths of the form ft = {l — t)g + tf where t G [0, 1]. It is a natural fact that the 
condition number plays an important role in the translation of the real-number arithmetic results 
to rational arithmetic results. We will see that this just produces an extra factor (the logarithm 
of the maximum of the condition number along the path) in the complexity bounds. Our main 
result thus reads as follows. 

Theorem 2 (Main). Assuming that g, f and zo are given in rational coordinates, and assuming the 
extra hypotheses \2.1 ) below, the homotopy method can he designed (see TrackSegment below) 



to produce an approximate zero with integer coordinates of f = f\ with associated exact zero 
Ci. The total number of steps is a small constant C^/nd^^'^ (d the maximum of the degrees of the 
polynomials in f , and n + 1 the number of variables) times Cq. The bit complexity of the algorithm 
is linear in Cq and polynomial in the following quantities: 

• n, S, d, h, where S is the number of nonzero monomials in the dense expansions of f and g 
and h is a bound for the bit length of the rational numbers appearing in the description of 

• log2(maxte[o,i]{/i(/t,Ci)}); o-^d 

• The quantity 



Here, \\ ■ \\ is Bombieri-Weyl norm (we recall its definition in Section 2.1) and fj,{ft, Ct) is condition 



number at the pair (ftXt), see (2.2) below. The bit size of the integer numbers in the coordinates 



of z^ is at most 0(log2(n) + log2/u(/, Ci))- 



The extra hypotheses that will be specified in (2.1) is, in words, that the angle between g 
and f — g is not too close to tt, that is that the segment {ft : t E [0,1]} does not point too 
straight-forward to the origin. 

The reader may note that, once this method is proved to work and programmed, it provides 
a status of mathematical proof to the path-following procedure. Moreover, its complexity does 
not differ much from the method of Theorem [Tj and the size of the integer numbers involved is 
controlled. 

In particular, the method can be used to give rigorous proofs to the results obtained via mon- 
odromy computation by algorithms in [28], [2]. Note that the implementations of the a-test carried 
out in ^20], [31] would not be sufficient for the above applications even in the situations where the 
zero count is known and the ct-test is capable of certifying the exact expected number of distinct 
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zeroes at the ends of the homotopy paths. This is due to a potential multiple path-jumping (which 
we just call path-swapping) resulting in a wrong permutation of zeroes produced by tracking mon- 
odromy loops. In other words, as we have already pointed out, while in certain cases the a-test for 
the end solutions can resolve the scenario of Figure [T| it is powerless in the scenario of Figure |2} 

1.1. Previous works and historical remarks. Smale [43] proved the first results on the average 
complexity of polynomial zero-finding using Newton's method. The problem of solving systems of 
polynomial equations became later one of the cornerstones of complexity theory in the BSS model. 
In the first paragraph of |13] we read "Also this work has the effect of helping bring the discrete 
mathematics of complexity theory of computer science closer to classical calculus and geometry". 
This paper is written with the intention of making another step in that direction. Theorem [2] 
establishes a strong link between the BSS model of computation and the classical Turing machine 
model, and we believe that this link strengthens both models at the same time. The Turing 
machine model is proved to accomplish a difficult task until now reserved to numerical solvers. 
An algorithm designed and analyzed in the BSS model is successfully translated (including its 
complexity analysis) to the classical model by carefully studying its condition number. Adding 
up, in the case of path-tracking methods for polynomial system solving, 

BSS model + condition number analysis = Turing machine model, 

the equality being strong in the sense that the complexity of the discrete algorithm is similar to 
the complexity of the BSS algorithm. This "translation" of computational models can probably 
be done for many of the algorithms originally designed in the BSS model. The reader may find 
many of our techniques useful for such a task. 

A natural precedent to our work is the algorithm in Malajovich's Ph.D. Theses [32] (see also 
[33], [M]) where a homotopy method for polynomial systems with coefficients in Z[i] is presented. 
The algorithm in [32] is certified under reasonable assumptions for e-machines (i.e. fioating point 
machines), in which some intermediate computations are made. Its total complexity is bounded 
by a number which does not explicitly depend on the condition number of the systems found in the 
path. In [32], [M] there are some "gap theorems" which give universal bounds for the value of the 
condition number in the case that the coefficients of the polynomials are integers (and assuming 
that the solutions are not singular). The biggest advantage of that approach is that a global 
complexity bound is found. By "global" we mean that it is valid for solving every generic system; 
as a drawback, that complexity bound is exponential on the bit size of the entries. 

In [HI Cor. 2.1], a universal upper bound for the bit size of rational approximate zeros of smooth 
zeros of systems with integer coefficients and degrees at most 2 is given. In general, one expects 
the condition number to be much smaller than in the worst case (see [IQ], [15] for results in this 
direction), and hence the bit size of the output of our algorithm will in general be much smaller 
than the upper bound of [HI Cor. 2.1]. 

One alternative to the rational arithmetic approach of this paper to the certification of homo- 
topies could be using interval arithmetic. For example, in a more general setting, [25] proposes 
step control by means of isolating a homotopy path in a box around an approximation of a point 
on the path. While the implementation of [25j does not provide certification, in principle, interval 
arithmetic can be used in an attempt to certify homotopy tracking using similar isolation ideas 
(see [21] for an ongoing work in this direction). 

1.2. Acknowledgments. In earlier stages of the ideas behind this work, we maintained many 
related conversations with Clement Fernet; thanks go to him for helpful discussions and comments. 
We also thank Gregorio Malajovich, Luis Miguel Pardo and Michael Shub for their questions and 
answers. Our beloved friend and colleague Jean Pierre Dedieu also inspired us in many occasions. 
The second author thanks Institut Mittag-Leffier for hosting him in the Spring semester of 2011. A 
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part of this work was done while we were participating in several workshops related to Foundations 
of Computational Mathematics in the Fields Institute. We thank this institution for its kind 
support. 

2. Technical background 

2.1. The vector space of polynomial systems. As mentioned above, we will center our at- 
tention on homogeneous systems of equations. For fixed n > 1 and an integer / > 1, let 
Til C C[XQ,...,Xn] be the vector space of degree / homogeneous polynomials with unknowns 
Xq, . . . ,Xn. As in IH], we consider Bombieri-Weyl's Hermitian product (■, ■)^; which preserves 
the orthogonality of different monomials and satisfies 

---^n )-^0 •••^n lUi- • 

For integers (ij > 1, 1 < i < n and an n-tuple (rf) = {d\^ . . . , (i„), we denote by 'H(rf) the vector space 
of systems of n homogeneous polynomials of respective degrees di, . . . , in unknowns Xq, . . . , X„. 
That is, 

The Hermitian product in ^{d) is then 

(/,fl') = {fl,9l)Hd^ H \- {fn,gn)Ha„ 

where / = (/i, . . . , /„), g = {gi, ■ ■ ■ , Qn) & T-i{d)- We also define 

11/11 = (/,/)^/'- 

The Hermitian product (norm) given by (■, ■) (|| ■ ||) in T-Li^d) is also called Bombieri-Weyl product 
(norm). It has a number of nice properties such as invariance under composition with a unitary 
change of coordinates, see for example [HI Section 12.1]. We denote 

§ = {/GH(,): 11/11 = 1}. 

Our main algorithm (Algorithm [l] below) will need to compute (■, ■) and || ■ p. This can obviously 
be done over the (complex) rationals Q[i], if the polynomials involved have coordinates in Q[i]. 

Remark 3. The extra hypothesis Theorem^ needs is 

(2.1) - Lo\\g\\ 11/ - g\\ < Re {g, f - g) < \\g\\ \\f - g\\, 

where Lq = 1 — 10^'^. While this hypothesis is not actually needed for the method to work (an output 
will equally be obtained if the hypothesis is not satisfied), it does simplify the complexity analysis 



of the algorithm significantly. Note that (2.1) means that the angle between the two vectors g and 
f ~ g in the space of polynomial systems is not too close to tt. This is satisfied by most reasonable 
choices of paths ft. 

2.2. Projective Newton method. We now describe the projective Newton method of [2S]- Let 
/ e n(d) and z e P(C"+i). Then, 

xp(/)(z) = z-(D/(^) ur'fiz), 

where Df{z) is the n x (n + 1) Jacobian matrix of / at z G P(C""*"^), and 

Dfiz) U 

is the restriction of the linear operator defined by Df{z) : C"+^ — )■ C" to the orthogonal comple- 
ment z-^ of z. The reader may check that Xp(/)(Az) = AXp(/)(z), namely Xp(/) is a well-defined 
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projective operator as long as the linear map Df{z) has an inverse. An equivalent expression, 
better suited for computations, is 



iVp(/)W 







where z* is the conjugate transpose of z. We denote by Np{fy{z) the result of / consecutive 
applications of Np{f) on initial point z. 

In general, one cannot expect that the solutions of systems of polynomials have rational coor- 
dinates. The goal of solvers is thus to produce rational points which are "close" in some sense to 
actual zeros. Following the approach of [H] , [H] we will consider that a point is an "approximate 
zero" of a system of equations if it is in the strong (quadratic) basin of attraction of the projective 
Newton method. Namely: 

Definition 4. We say that z G P(C"+^) is an approximate zero of f E 'H^d) with associated (exact) 
zero C G P(C"+i) if Nv{f)\z) is defined for all I > and 



d^(iVp(/)'(.),C)<^^^ 



22 

Here dji is the Riemann distance in P(C"+^), namely 



/ > 0. 



\(z z')\ 

dniz.z') = arccos^pi^ G [0,7r/2], 

where (■, ■) and || ■ || are the usual Hermitian product and norm in C"'^"'^. Note that dii{z, z') is the 
length of the shortest curve with extremes z, z' G P(C"'^^), when P(C"'^^) is endowed with the 
usual Hermitian structure (see for example [HI Page 226].) 

2.3. The condition number. The condition number at {f,z) G 'H(d) x P(C"'''^) introduced in 
jH] is defined as follows: 

(2.2) M/,^) = 11/11 \\{Df{z) \,.y'Biag{\\zt-'dy')l 

or n{f,z) = oo if Df{z) | ^± is not invertible. Here, ||/|| is the Bombieri-Weyl norm of / and the 
second norm in the product is the operator norm of that linear operator. Note that /i(/, z) is, up 
to some normalizing factors, essentially equal to the operator norm of the inverse of the Jacobian 
Df{z), restricted to the orthogonal complement of z. Sometimes fi is denoted /inorm or /iproj, but 
we keep the simplest notation here. One of the main properties of /i i^that it bounds the norm 
of the implicit function of the mapping (/, z) (-)■ f{z). In other words, following [11], Sec. 12.3 and 
12.4]: 

Lemma 5. Let {g, Co) £ 'H(d) x P(C"+-'^) be such that g{Co) = 0, and fi{g, Co) < oo. Let ft-, t G [0, e) 
he a curve in l-Li^d), fo = 9- Then, for sufficiently small t < e, (q can be continued to a zero (t of 
ft, that is there exists a curve t h-)- C(t) C P(C"'+^) such that C(0) = Co o.^'^d, denoting ({t) = Q, 
we have ft{Ct) = for every sufficiently small t. Moreover, the tangent vectors satisfy: 

llColl <M^7,Co)||/o||. 



We will also use a variation of this condition number, namely Xi in equation (3.4) below. 
The following result is a version of Smale's 7-theorem (cf. [S]), and follows from the study of 
the condition number in 



""^This property inspired its definition, see [H]. 
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Proposition 6. [H Lemma 6] Let C e P(C"+^) be a zero of f E 'H(d) and let z G P(C"+^) he such 
that 

dniz, C) < ,0/0 % >.N , where Uq = 0.17586. 
Then z is an approximate zero of f with associated zero (. 

2.4. Complexity and the condition metric. According to [31], the complexity (dominated by 
the number of Newton steps or number of while loops) of an algorithm performing the homotopy 
method should depend on the so-called condition length of the homotopy path. Given a path 
(htXt), to <t <ti where (t is a zero of ht and ht G P('H(d)), (t € P(C"+^), the length of the path 
{ht, Ct) in P(H(d)) X P(C"+i) is given by the integral 

dt 

'r('^t,Ct)(P(«(d))xP(C"+i)) 







/" 

J to 





to ' 

One must here understand ht and (t as tangent vectors in Thtf{T-L{d)) and Tc^^FiVJ^^^) respectively. 
If ht and C,t are given in coordinates, this means: 

ii;„2 _\\htr \{Kht)? 

Wilt 
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Now, the condition length (or length in the condition metric) of the same path is defined in |39] 

as 

/^(^i' Ci) Y ll^t|lThjP(-H,^)) + IICt|lT(;jP(C"+i) dt. 

to 



Note that, from (2.2), fi{htXt) only depends on the projective classes of ht and Ct and this last 
integral is thus well defined. 

Now, given a path {ft, (t) £ 'H(d) x P(C""'"^) where ft{Ct) = 0, we define its condition length as 




(2.3) Co = Co(/t,Ct) = 

That is, ColftXt) is the length in the condition metric on S x P(C"^^) of the path obtained by 
projecting {ft,Ct) £ ^(d) x P(C""'""'^) on S x P(C"'"'"^). The reason of this change is that segments 
in 'H(d) project nicely on the sphere S (indeed, they project onto pieces of great circles in S), but 
they do not project so nicely on P('H(d)). This makes our analysis easier and, presumably, has 
little effect in the complexity bounds. 

Note that, if ft G H(d) is the horizontal lift of ht G P(7{(d)) then the condition lengths of {ft, (t) 
and {ht, (t) coincide. In the general case, however, the condition length of {ft, Ct) is greater than 

2 

that of (7rp(^^^j)(/j), Ct), because in general Re {{ft, ft)) < |(/t, ft)\^- 

3. A ROBUST HOMOTOPY STEP 

In this section we set up the backbone of our main algorithm: how to correctly choose a homotopy 
step. We do this by stating a theorem that, given sufficiently close polynomial systems g, f E "H^d) 
and an approximate zero zq of g associated to an actual zero Co of 5'; guarantees that Co can be 
continued to a zero Ci of /, and moreover a projective point sufficiently close (in a sense that we 
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will precisely determine) to Np{f){zQ) is an approximate zero of / with associated zero Ci- There 
are some precedents to this result in [32] but our theorem is needed to get the sharp complexity 
bound of 



3.1. Some constants. As it is common in the exphcit description of many numerical analysis 
algorithms, we will need to use some constants, that need to be described explicitly because they 
intervene in the definition of the algorithm. We will use a free parameter 1/2 < 5 < 1 that will be 
set to 3/4 in our implementation of the algorithm. The rest of the constants are: 

0.17586 (the constant from Proposition |6]) , 



P = V2+ + 

(3.1) c' = c's 

Finally, let c be any number satisfying 



{25 - l)uo 1 

a = as = —p^ < —;=, 

V2 + 25uq ^2 



a)^ < 1. 



(3.2) 



C < C5 



c . 



1 + V2uq/2 

Only the value of c^/(2P^) will appear in the description of our main algorithm (see Algorithm 
TrackSegment below.) Hence, one can choose any value of (c^/2P^) G Q such that 



0.00034412... 



2P2 - 2P2 

In our algorithm, we will choose iMf^ = 0.00034. The reader may check that the following holds. 



(3.3) 



Pfl 



< 



d 3a 
P - 2V^' 



3.2. A version of the condition number. The condition number jJ^ifX) of (2.2) above can be 
computed using a more amenable expression if f{() = 0. 

Let g,g & 7i{d) be two polynomial systems and let z G P(C"'"'"^). Let Xi = Xiid^^): Xi = 
X2i9, 9, z) and (f = (p{g, g, z) be defined by 



(3.4) 



Dg{z] 

z* 



(3.5) 



X2 



[\fdi\\g\ 



\ 



Ii2 \\9\ 



Dg{z] 

z* 



-1 



9\z, 




1/2 



(3.6) 9? = x\X2- 

Note that these formulas do not depend on the representative of z and thus are well defined. Their 
value is also invariant under multiplication of by a non-zero complex number A G C. 

It was noted in % eq. (2.2)] that if t ^ {ft, Q) CSx P(C"+^) is a curve such that ft{(t) = 0, 
then 

(3.7) xiiftXt) = KftXt), viftJtXt) = KftXtMftXt)]], vt. 

Thus, Xi(/, z) is a version of the condition number fj,{f, z) (equal if f{z) = 0) and ip is, if {ft, Ct) ^ 
§ X P(C"+^) and ft{Ct) = 0, the quantity inside of the integral defining Cq in (2.3). 
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From Lemma [sj with the notation of ( 3.7[ ) we have: 



(3.8) ¥'(/o,/o,Co) <M/o,Co)||/o|IVl+M/,Co)2- 

The reason to use Xi instead of ^ and (p instead of simply the term inside the integral defining Cq 



is that the rates of change of Xi &nd (p are easy to analyze, see lemmas 17 and 18 in Section 10 



3.3. A robust homotopy step. We now state our main technical tool, which is a more detailed 
and complete version of Lemma |5] about continuation of zeros, designed to answer the following 
questions 

• how long can a zero of g be continued when g is moved? 

• how do Xi cLnd (p vary in this process? 

• if an approximate zero of g is given, how long will it still be an approximate zero as g is 
moved? 



We will use the constants defined in Section 3.1 Given two systems f,g ^ T^{d), we consider the 
Riemannian distance in the sphere § from (?/ to //||/||, that is: 

, , 9 f \ ^e{g,f) 
as 7^— 77, 77-777 = arccos -— 



^11 11/11/ ll^l 

Theorem 7. Let g, f ^ '^{d) be two systems of polynomial equations such that g ^ Xf "iX E 
Let zq be an approximate zero of g satisfying 

for some exact zero Co of g. Let 

\\grf~Re{{f,g))g 



9 



MV\\m9r-M{f,9)r 

That is, g is the derivative at of the arc-length parametrized short portion of the great circle in 
S, fromg/\\g\\ tof/\\f\\. Let 

Xi = Xi{9,Zo), X2 = X2{9,9,Zo), (f = ip{g, g, zq). 

Assume that 

(3-10) ds(A^Al)< ' 



grwfwj - Pd^i^ip- 

Then, 

(1) Co can he continued following the straight line homotopy 

(3.11) f, = {l-t)g + tf 

to a zero Q of ft, namely, there exists a curve t ^ ({t) such that ({0) = (q and, denoting 
at) = Ct, we have ft{Ct) = for t e [0, 1]. 

(2) We have the following inequality: 



(3.12) ^ < 



V2fi{g,Ci 







\2 



;i - V2uo/2) 



1+V2'' 



(3) The condition length Co{ft,Ct) of the path (ftXt) defined in (2.3) is essentially equal to 
^ds yyj ■ More exactly: 

p (1- V2V2)'+-/^ln(l + e') ^ Co(/.,C.) ^ 1 + 
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(4) For every z E P(C"+^) such that 

(1 - S)uo 



(3.14) rf«(z,iVp(/)(zo)) < 

we have that 

(3.15) dR{2, Ci)< 



we have that 

Uo 



2rf3/2^(/,Ci)- 

In particular, z is an approximate zero of f with associated zero (i . 



The proof of this theorem is a long and tedious computation. We delay it till Section [10 

4. A SCHEMATIC DESCRIPTION OF THE ROBUST LINEAR HOMOTOPY METHOD 

In this section we describe an algorithmic scheme for the linear homotopy method. The procedure 
in this section is not quite an algorithm, because we do not specify how to perform some of the 
tasks it requires. We will however prove that any actual algorithm designed to fit into the scheme 
of this section has certified output and the number of iterations it performs is essentially bounded 
by the condition length Cq. 

Let g, f be two non-collinear systems, that is g ^ \f for every A G M. Let ft be defined by 
( |3.1l| ), so /o = g, fi = f- Let S eQ, R E ^M with 1/2 < 5 < 1 and > ^2 be two arbitrary 



constants B 

Algorithmic Scheme 1. = TrackSegment_Scheme(/, gf, zq) 
Require: f,gE 'H{d) non-collinear with coefficients in Q[i] 



zq E Q[i]'^~^^ is an approximate zero of g satisfying (3.9). 



Ensure: z^, E is an approximate zero of f associated to the end of the homotopy path 



starting at the zero of g associated to Zq and defined by the homotopy (3.11). 
1: i ^ 0; Si ^ 0. 
while Si ^ 1 do 

9i ^ fs, ■ 
Let 

\\giff-Re{{f, gi))gi 



9i 



\gi\W\\f\ngi\? -Re{{f,giW 



LetXi,i = XiA9i,Zi), Xi,2 = XiAdi, 9i, Zi) and (fi = ipi{gi, gi, Zi) as defined in (3.4), (3.5) and 



(3.6). 



5: Let ti be any positive number such that 

4.1 L < " ^^ '■^ < Ur 

- \\g.\W\M' + 2t. Re {g^, f - g) + mf - ~ 



where 



_ 



If such ti does not exist, let 1 (which will imply that the while loop finishes in this step). 

6: if ti > 1 — Si then 

7: tii-l-Si. 



^If R > \/2 then values ti as described in Algorithm [l] exist. They also exist for smaller values of i? > 1 like 
R = 1.0003 but taking i? > ^2 will make our formulas look prettier. This assumption does not affect much to the 
running time of the algorithm. We will only use R^ in the algorithm. Hence, we can take R G 
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8: end if 

9: Sj+i ^ Si + ti; 

(1 - 5ful 



10: e i- 



11: 



12: ^ any vector in Q[i]'"+-'^ satisfying 

(4.3) dR{zi+i, ZiJ^i) < ^/e 



[l - 5)uo 



2c/3/2(l + 35V2)Xi,i' 



13. 
H 

15 
16. 



i ^ i + 1. 
end while 

2* ^ ^i+i- 



Note that (2.1) implies that for i > we have 
(4.4) - Lolkill 11/ - QiW < Re {gij- Qi) < \\gi\\ \\f - gi\\, 

where Lq = 1 — 10~^. 

Theorem 8. The output of any algorithm performing the instructions described in TrackSeg- 
MENT_SCHEME is certified. Namely, for every i > 0, the point Zi is an approximate zero of 
gi = fsi, with associated zero^C,i, the unique zero of gi such that {gi,(i) ^^^s in the lifted path 
iftXt)- Moreover, 



Let Cq be defined by (2.3) and (3.11). If Cq < oo, there exists k > such that f = g^. For the 
number of homotopy steps k the following bounds hold: 

C'd'/^Co<k< \Cd^/^Co], 

where 

(J ^ ^,^P 

c(l-y2uo/2)i+^ln(l + c')' C 
In particular, if Cq < oo, there exists a unique lift {ft,Ct) of the path ft, and the algorithm finishes 
and outputs z^,, an approximate zero of f = g^ with associated zero (k, the unique zero of f such 
that (/, Cfe) li^s in the lifted path {ft,Ct)- 

Finally, the two following inequalities hold at every step of the algorithm: 

V2fi{g.,Qf 



(4.5) ifi < 



;i - V2uo/2) 



1+V2'' 



(4.6) Ur — L > — — — — , where c is a universal constant. 

d-^ maxt6[o,i]{/x(/t,G)*} 

Remark 9. As said above, we will choose 6 = 3/ A in our main algorithm TrackSegment. With 
that choice and the use of Frobenius norm instead of operator norm for the computation ofxi,i, 
our practical implementation we will have 

28<C' <C < 79Vn + 1. 
^Note the slight abuse of notation: we just use Q the zero of gi — fs-, so we should actually denote it by Cs;. 
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Note that C is not a universal constant as it depends on n. The value of c is needed only for 
the bit-complexity analysis where it will be replaced by an 0(1). One can however estimate it as 
c ^ 0.00003. 

Proof. The proof of correctness of the algorithm is by induction on i. The base case of our induction 
i = follows. Assume that 

'"'^'"'^'^ - Mvfc)- 
We claim that we are under the hypotheses of Theorem [7j Indeed, 

9t+i = fs,+i = fs,+u = (1 - Si - ti) g + {si + ti) f = 
(1 - s^)g + Sif - tig + tif = gi + ti{f - g). 

Thus, 

Re{gi,gi+i) ^ WgiW^ + tjRe {gjj - g) ^ ^ 

UiW ^/M^ + 2tiRe {giJ - g) + t^\\f - g\\' ~ 
Thus, from Lemma 10 below we get 

rAQ\ A ( 9i 9i+i \ Re (5fi,5fi+i) c 

(4.8) ds - — = arccos - — — < arccosL < — . 

WmW ll5'»IMl5'i+i|| Pd-y^'^i 

In particular. Theor em [7| applies to the segment proving our induction step and also 

proving (4.5) from (3.12). Additionally, if the i-th step is not the final step in our algorithm 
(equivalently, (/j+i 7^ / or s^+i < 1) then we have 

Re {gi, gi+i) \\gi f + ti Re {gi, f - g) 



\\gi\\^M^+'^kRH9U^W+nW^^ 

which again using Lemma [TO] yields: 



< U. 



fAn\ ^ I 3i 9i+i \ ^^{gt^gi+i) ^ c 
(4.9) a§ I — f7, = arccos - — — > arccos (7/? > 



9if\\9i+i\\J \\9i\\ \\9i+i\\ RPd^l'^Lpi' 

Now we prove the bound on the number of steps. From (3.13) and (4.9) we have that, as long 
as Si+i < 1, 

V IL^lF ILAF ik<IItc,p{c"+i)"^ - 

( g, g,+, \ (1- v/2no/2)i+^^ln(l + cO 

'fids ]]— f7, 71 [7 > 

Vll^ill \\9i+i\\J c 

c(l- v/2txo/2)i+^ln(l + cO 
c'RPd^/^ ■ 

Thus, as long as Sj+i < 1, we have 



WftW Re ((/,/*))' 



Co - I Jf^ |[^jj4 ^ IICJt^^p(C"+1)^^ > 



V r\(f nt/l'^*"' UP dt> 
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{i + l)c(l - V2uo/2y+^ln{l + c') 

In particular, if Sj+i < 1 then 

i + 1 < ^ P Co. 

c{l-V2uo/2y+^ln{l + c') 

The first non-negative integer i which violates this inequality is thus an upper bound for the 
number of iterations of the algorithm. This finishes the proof of the upper bound on the number 
of steps. For the lower bound, note that, even if Sj+i = 1, from (3.13) and (4.8) we have 



""uff nJI'^*"' ^^iitM! , ||.||2 dt< 

gi Qi+i \ l + ^uo/2 c(l + y2uo/2) c' 
fids jr ii V. ,7, ,^^./o < TTI^TTT; ,^..r. < 



kill' lk*+illy (1 - V2no/2)^ Prf3/2(i _ y2?zo/2)^ " Pd^''"' 
Thus, if k is the number of iterations needed by the algorithm (i.e. Sfc_i < = 1) then 



1 /llf.l|2 " "-^2 



r\(f cu^^^'^^' dt<k 



i=0 *J 

In particular, we conclude that the total number of iterations is 

k > , 

c' ' 

which is the lower bound on k claimed in the theorem. 



For (4.6), note that 

Ur-L= .r^::^^ i 1 



Using R > \/2 and roughly bounding the term inside the parenthesis, we get 



Ur-L> _ J',^ . > 



,2 c2 f(l-y2W2)'+^^' 



which implies (4.6). 

Lemma 10. Let L, Ur he defined as in the algorithm. Then, 

> Rph^y ™^ < p¥r^- 

Proof. We prove the second inequality. Recall the elementary fact that for < s < 1 we have 

cos(s) < 1 1 . 

^ ' 2 24 

Then, because arccos is a decreasing function in [0, 1], 

(4.10) arccos ( l-- + —j <s, se(0,l). 



□ 
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In particular, 



c2 



arccos L = arccos 1 — ^ ,„ „ + ^ , ^ . ,„ . < 



as desired. The first inequality is proved in the same way, using that for < s < 1 we have 

cos(s) > 1 - 



5. Computational considerations 



□ 



Provided Theorem [8| the rigorously certified homotopy tracking could be accomplished by 
way of exact rational arithmetic employed in all of the computations described in TrackSeg- 
MENT_SCHEME. In this section we discuss some of the aspects of this issue, to facilitate the reading 
of our main algorithm TrackSegment below. 

5.1. Operator norm vs. Probenius norm. In Step 4 of TrackSegment_Scheme we need 
to compute the operator norm of a matrix, which is a non-trivial task. Actually, one just needs 
to the square of such norm, to use it in steps 5 and 10. Instead of computing the square of the 
operator norm, one can just compute the square of the Frobenius norm IKajj)]!!' = 
which involves only rational computations. Both norms are related by the inequalities 

ll-f <l|-|||<(^ + l)l|-f- 

On the other hand, which is the squared norm of a vector involves only rational computations 
and thus can be computed exactly. Then, instead of iff in Step 5 we can use the product of 
and a version of using the squared Frobenius norm. 

Let us put this in a general framework. Assume that we can compute some quantity xfi 
satisfying ^ < x\i < S'^xl,! for some S > 1. Let L, be computed with the same formulas as 
L, but using Xj i instead of Then, it is easy to see that L > L and 

~ _ _ c2 c2 _ 

~ ^ ~ - ^ " AS^PM^^f ~ 

Thus, if we find ti such that 

hif + {QiJ - g) 

- MVUW' + 2t. Re {g., f-g)+tnf- g\\' - 
then in particular the hypotheses of Theorem [s] are fulfilled changing R to \/2S and the number 



of steps is at most 



k < \ ^ ^^^t d^/'Co]. 

" c(l-y2uo/2)i+^ln(l + c') ' 



In particular, we have proved the following. 

Lemma 11. If in TrackSegment_Scheme, R is changed to \f2 and Xii is changed to 
(defined the same way as xli but using Frobenius norm instead of the operator norm), then any 
algorithm performing the computations in TrackSegment_Scheme has certified output in the 
sense of Theorem\^ The number of iterations is at most 



k<\ Vnn + l)Pc' ^3/2 ^ \79V^d'/'Co]. 

~ c(l - V2uo/2y+^ ln(l + c') ^ouh 5=3/4 ' ' 
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Moreover, at every step we have 

IC)--LJ Xi,lXi,2 S 



;i - 72^0/2)1+^ ' 

^^•^^ ^V^(^) nd? max{/i(/i, CtY} ' 

where c is a universal constant. 

5.2. Computing the step size. Step 5 of TrackSegment_Scheme requires finding a t sat- 
isfying (4.1), wliicli a priori means computing approximately the smallest positive roots of two 
quadratic polynomials. An easier way to get this is using bisection method to locate a root of the 
equation 

Oi + te2 L + u ^ 

a(t) = — = = 0, 

y/ei{e, + 2te2 + t^es) 2 

with stopping criterion given by 

, / X, U-L 

Ht)\ < ^— . 

If t G M satisfies this stopping criterion, then (taking 9i = WgiW^, O2 = Re {gi, / — (?), 6*3 = ||/ — gW^) 
it also satisfies (4.1). When applying bisection, we need to be able to determine the sign of a{t). 
It is not hard to accomplish this without computing square roots using the following subroutine. 

Algorithm 1. (s,r) = Computesign(6'i, 612, 613, t, L, f/) 

Require: 9i, 92,63, t, L,U E Q, 61,63 > 0, 0| < 6163, < L < U < 1. 
Ensure: s = 1 if a{t) > 0, s = — 1 otherwise, and 

{6i + t62y 

e,{6, + 2t62 + tW3) 

{6i + t62? 

• ^ 6^{6i + 2t62 + t^63) 

2: if 61 +t62>0 and r > (L + Uf/A then 
3: S ^ 1 

^; else 

5: Si 1. 

6: end if 

The bisection method mentioned above is then as follows. The requirement 1 — 10^^ < L in 
the description of the following algorithm corresponds to the extra hypotheses (2.1) in our main 
algorithm. 

Algorithm 2. t = LUquadratic(^i, ^2, ^3, -f', f^) 

Require: 61,62,63, L,U e Q, 61,63 > 0, -1^/6^3 < 62 < ^fWz, 1 - 10"^ <L<U <\. 
Ensure: t = f e Qn [0,1], m e Z, < I e Z such that 

(5.3) L < , + = < U, 

- y/e,{6, + 2t62 + tW3) - 

if such t exists (otherwise, output t = 1), and such that 

(5.4) < m < 2' < max 1 



6x{lJ-L) 
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10. 

11 

12. 

13: 

14: 

15: 

16: 

11: 

18: 



(si, ri) ^ Computesign(6'i, 02, 0^, ti, L, U) 
if 01 + 02 > and ri > L2 then 

t ^ ti 
else 

U2 ^ 

to + h 

to i — 

2 

(52,^2) COMPUTESIGN(6'i,6'2,6'3,t2,^, 

/ ^ 

while L2 > r2 or U2 < r2 or 0i + t202 < do 
if S2 = 1 then 

to ^ t2 
else 

tl i- t2 

end if 

to + tl 
to < z 



19: {s2,r2) ^ Computesign(0i, 02, ^3, ^2, L, U) 

20: I ^ l + l 

21: end while 

22: t ^t2 

23: end if 



Lemma 12. LUquadratic produces t = m/2' satisfying (5.3) and (5.4), or t = 1 in case there 
exists no t that satisfies (5.3). Moreover, the number of iterations it performs is at most 



O l0g2 



0^ 



Proof. Let 



We first claim that 



(5' it) = a\t) 



0i{U-L) 

01 + 102 



^0i{0i + 2t02 + t^0^ 



{01 03 - 02') t 



^l{03t^ + 2 02t + 0i 



This is a routine computation and is left to the reader. In particular, 6*1 < 0i0^ implies that a{t) 
and /3{t) are decreasing functions for t > 0. If > L (which is decided in Step|4]of Algorithm]!]) 
then there are two possible scenarios: 



(1) If /3(1) > U then a t satisfying (5.3) does not exist and the output of the algorithm is t = 1 
as claimed. 

(2) If /3(1) < U then the output of the algorithm t = 1 satisfies (5.3) and (5.4) as claimed. 
On the other hand, if /3(1) < L, then 



L + U 
a{0) = 1 — > 0, 



a(l) = /3(l)- 



U L-U 
— < — - — < 0, 
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which imphes that the bisection method used in the algorithm produces an approximation of the 
unique root G (0, 1) of a(t). In particular, note that = implies that 9i + 1^92 > and 

L + U 



which yields f3{t^:)^ G (L^,f/^). By continuity of (3, we conclude that the algorithm will at some 
point compute a t2 such that L"^ < r2 < U"^ and Oi + ^2^2 > 0. That is, the algorithm finishes at 
some point, and the output satisfies (5.3) as claimed. It is a simple induction exercise to prove 
that, whenever the condition of Line 12 is satisfied (that is to say, at every step of the algorithm, 
except possibly at the last one) we have 

b,g] c [to,ti], 

where 

[p,g] = {tG[0,l]:|a(t)|<(f/-L)/2} 

= {tG[0,l]:/3(t)G[L,f/]} = r'([^,f/]). 
At every step of the algorithm, the bisection method satisfies 

1 

Thus, if the /-th step is not the last step of the algorithm then we have 

1 



q — p < 



2'' 



From the Mean Value Theorem of calculus, we have 

U-L 



q-p 

for some i G [p,q] = (3~^{[L,U]). Thus, 

1 



\q-p\ 



q-p 



(5.5) 2' < 

Note moreover that 
(5.6) 

Now we have to distinguish two cases: 



P'{t) ^max{|/3'(t)|:/3(t)G[L,f/]} 



U-L 



U-L 



h03 - 0^ 



(1) If ^2 > then and t < 1 yield 

(2) If ^^2 < then by hypotheses we have 62 = — e with < e < Ly/6^. Thus, 



2m ^ ^3 



/3 



ei/2 



ly 4e2 



and hence /3(t) > L implies that t < 6i/{2e). Thus, (5.6) implies that 
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which readily gives 



\m\<P%tff<f, te[p,g]. 



Thus, for every possible value of 62 G {—Ly/OiOs, we have that 

max{|/3'(t)|:/3(t)G[L,f/]}<^. 

C71 



This together with (5.5) proves that, if the l-th step is not the last step of the algorithm, we have 

2' < 



»,([/- L)' 

For the last step, this quantity has to be multiplied by 2. The last claim of the lemma follows. □ 

Remark 13. Our implementation of LUquadratic continues bisection if the denominator of its 
output t is larger than the denominator of Si in step 23 0/ TrackSegment (Algorithm^ until 
the denominators match. This is done in order to reduce the size of the denominator o/sj+i. 



5.3. Finding a close-by number with small integer coordinates. In Step 12 of TrackSeg- 
ment .Scheme we change Zj+i to a close-by vector Zj+i with rational coordinates. Although Zj+i 
already has rational coordinates, we need to replace Zj+i with a nearby vector whose coordinates 
are integer numbers of bounded (small) absolute value. If this step is not performed, the number 
of bits required to write up ZiJ^i might increase at each loop, which is to be avoided. In this section 
we show how to deal with the general problem of, given z G Q[i]"^^ and e G Q, e > 0, finding 
z G Z[i]"~''^ such that|^ z) < ^/e, and such a way that the absolute value of the coordinates of 
z is relatively small. 

Let us consider the following algorithm. 

Algorithm Z. z = ShortZero(2;, e) 
Require: z G Q[i]"+^■ e G (O, \) fl Q. 
Ensure: z G TJ^iY^^ such that 

(5.7) dR{z,z)<y/e, 

and such that the integer numbers appearing in the expression of z are bounded in absolute 
value by 3i ^ 



e 



1: Let — + i-^, ai,bi,Ci,ei G Z, Cj, > 0, < i < be the coordinates of z. 
2: {co---Cn) ■ (eo---e„). 

3: X ^ m ■ z. 



4: r ^ 



2V ^ 
20 



9. 
10. 
11 



a ^ A 

while a < do 

k ^ k + 1 
end while 

z ^ [2"^x] . 



^Recall that dji{x,y) = arccos n^'^i'if^i is the usual distance from x to y as projective points in iP''''r"+i 
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Here, by [y] {y G C"+^) we mean the following: if y = (oq + i&o; • • • , c-n + i^n) then 

[y] = (M + i[bo], . . . , M + i[bn]), 

where for t G M, [t] is the integer number which is closest to t and is smaller than t in absolute 
value (that is, [t] = [tj if t > and [t] = \t] if t < 0). 

Lemma 14. Let < 63 < 6i. Then, the function 

has a global minimum value equal to a/1 — ^s/^^i- 
Proof. Note that w is a differentiable function and 

^3 + ^2 



Hence, the minimum of w is attained at 62 = -\f(hOl, 02 = \f(hh or 62 = -63- Now, ^(^^3) = 
w{-VM^) = 1 and ^(-^3) = Vl - 03/61 < 1. The lemma follows. □ 

Lemma 15. Algorithm^ produces z = (ao + if3o, . . . , a;„ + i/3n) G Z[i]"+^ satisfying (5.1) and such 
that 

\a^\,m<3J^- VO<z<n. 

Proof. First note that, if the stopping condition of the loop is satisfied at the first step, then the 
output z = X of the algorithm satisfies dji{z, z) = and 



e \ e 

and hence the claim of the lemma follows. Otherwise, the numbers a, k computed by the algorithm 
satisfy 

II II 2 

^ ^ ' - 2{n + l)r 

Let X = (ao + i/^o, • • • , a™ + i/3n) be the coordinates of x. Then, for z = 0, . . . , tt,, we have: 

|(«, + iA) - 2-\a, + = (p-'^a,] - 2-'a,)^ + (p-'^A] - 2-^=^)2 < 2. 
Hence, denoting y = 2~^x and v = z — y we have 

n 

(5.9) \\vf = 5^ |(«. + iA) - 2-^«. + iA)|' < 2{n + 1). 

i=0 

On the other hand, 

|2 ||„.||2 11 .||2 ||~ „.||2 r,o„/;: „.\ II=I|2 



\vr = \\yr~\\z-yr = 2Re{z,y)-\\z\ 
2 (^2-'=a,[2-%]+2-^'A[2-'=A] 



. j=0 



>2 K^[2-'^a,]2 + [2-'=/3,]^ 

\i=0 

= 2||5|p- pf > 0. 



\zf 
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That is, WvW^ < Hence, the use of Lemma 14 in the following chain of inequalities is justified: 

\{y + v,y)\ ^ ^e{y + v,y) ^ \\y\\^ + Re {v,y) ^ 



\\y\\^ + Re{v,y) 



\y\\\\y + '"\ 



Thus, 



> 

I/||\/|bP + 2Re {V,y) + ||f P LemmaM 

dniz, z) = dii{z, x) = dniz, y) = dniy + v,y) = arccos 1/^, j" ^'^11 < 

\\y\\\\y + v\\ 

2''\\v\\ 2^J2{n + l] 

< arcsin — 

5:1 



arccos W 1 — 



arcsm 



arcsm 



Note from ([5^ that 
(5.10) 



X y r 2 



The reader can check that the function s h> s ^ arcsin s, s G [0, 1) is an increasing function. From 
this fact and (5.10) we get: 

2*^2(n+l) 



arcsm 



2'=^2(n+l) 



< 



arcsm 



< 



21 
20 



which readily implies 



. 2V2(r2 + i) ^ 2y2(!iTT; 

aij(z, < arcsm — ^^-7^ — < vr — — 



\x\ 



\x\ 



01 



as wanted. 

For the bound on I a,- 1 note that 



[2-'''ai] \ < \2-'a,\ < 



A/8(n + l)r 
Ikllv^ 



aA < 3 



n 



where we have used ||a;|| > \ai\. An identical chain of inequalities works for /3j. 



□ 



6. The main algorithm 

We now describe the pseudo-code of an actual algorithm that performs the instructions described 
in TrackSegment_Scheme and is thus certified. 

There are two choices in TrackSegment.Scheme: R and S. We choose R = \/2 and 5 = 3/4 
which make the computations simple. Besides, instead of usi ng o perator norm for the computation 
of xli we use Frobenius norm, which according to Section 5.1 multiplies by a factor of \/n + 1 
the upper bound for the number of homotopy steps. The reader may find helpful Table [2] for 
comparing the names of the variables in TrackSegment.Scheme and TrackSegment: 

Algorithm 4. 2^, = TrackSegment(/, ^f, 2:0) 



in+l 



is an approximate zero of g satisfying (3.9). 



Require: f,g E H^d); e 

Ensure: z^, G Z[i]"~^^ is an approximate zero of f associated to the end of the homotopy path 
starting at the zero of g associated to Zq and defined by the homotopy (3.11). 
1: i ^ 0; Si = 0. 
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Table 2. Notation of TrackSegment and TrackSegment.Scheme. 



Tr agkSegment 


TrackSegment Sgheme 


III 


II fp 

\\J II 


112 


lloF 

Hi/ II 


71 o 


Re (f a) 
^'^'^ \j 1 y/ 


fl 


II f - olP 


T14. 


llolP 
11^*11 


715 


Re (/, Qi) 


ne 


Re (/ - g, Qi) 


ny 




Vl 




V2 


gi{.Zi) 


M 


mgi{zi)\-^ 

\ Zi / 


M 




a 


xi {— xi computed with Probenius norm) 


b 


xi 


ab 


= Xixi (plays the role of (p'^ = xixi) 


W 





2: n, ^ Wff. 

3: na ^ \\g\\^. 

4: na ^ Re{f,g). 

5: h nl + nl — 2n^. 

ul 

^ (4d)-^(l + 9V8)' 

while Si < 1 do 
n4 ^ (1 - Si)^n2 + sfrii + 2sj(l - Si)n3 
n5 (1 - Sj)n3 + s^ni 
ne Sjni - (1 - Sj)n2 + (1 - 2sj)n3 

^ Ikip 



5, 

10. 

11 

12. 
13. 
H 



15: 



16. 
17. 
18. 
19. 



Ml ^ Dg{zi); ^ Df{zi 



/ X /(I - Si)Mi + SiM2\ 



n n— 1 



k=0 1=0 / \fe=0 



^2 ^ 9i{zi) = (1 - Si)g{zi) + Si^i e 
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20: 



21 

22 
23. 
24 
25. 



W ^ Wo/{ab) 

L^l-W + Wy6; U ^1- W/2. 
ti LUquadratic{n4^, uq, h, L, U) . 
Si+i ^ min{l, Si + tj}; 
e ^ So/a 



26: M^(^^-^^+^)^: + ^^-^^^^^)eA^„^,(C). 



27. 

28. 
29. 
30. 
31 
32 
33. 



V5 ^ 9i+iizi) = (1 - Si+i)g{zi) + Si+iVi e C 

Zi+i ^ SHORTZERO(2;j+i,£:). 

i ^ i + 1. 
end while 



We should point out that in our practical implementation of the algorithm lines 13, 14, and 
19 as well as lines 26, 27, and 28 correspond to the calls to the subroutine executing one step of 



Newton's method for a specialization of the system (3.11). We break this up into smaller steps 



above for the purpose of the complexity analysis performed in Subsection 7.4 



Remark 16. From (5.2) and (5.4), for every i>0 the number of iterations o/ LUquadratic at 
Step 23 is at most 

||/-(7||nt/3max{/i(/i,Ci):0<t<l}~ 



O loggmax 1, ■ n\f\\ a ^ + ^ ^\ 

V V mm{||/t|| : < t < 1} 

7. Complexity analysis 

In this section we analyze the bit complexity of TrackSegment. Given a rational number 
p/g G Q, gcd{p, q) = 1, the bit length of p/q is defined as 

h\{p/q) = log2(max \p\, \q\) + 1. 

We also define bl(0) = 1. Note that bl(p/g) is a (tight) upper bound for the number of binary 
digits required to write up p or q. writing p/q thus takes at most 2bl(p/q) bits. 

Recall that an algorithm (i.e. a Turing machine) is said to have running time polynomial 
on quantities ci(x), C2(x), . . . , c/(x) (where the q(x) are quantities depending on the input x of 
the machine) if there exists a polynomial p{X) G M[Xi, . . . , Xi] such that the running time of the 
machine on input x is bounded above by p{ci{x), . . . , q(x)). A convenient notation is the following: 
given some function f{x) depending on the input x, we say that 

/(x)<(ci(x),...,q(x))«« 

if a polynomial p exists such that f{x) < p{ci{x), . . . , q(x)) for all possible input x. If a machine 
has running time which is polynomial in the (bit) size of its input, that is if the running time of 
the machine is input ^size'^^^^ then we say that the machine works in polynomial time. The reader 
does not need be very familiar with the concepts of computational complexity or Turing machine 
model to understand this section. However, we quote [TTl, Introduction] and its references for a 
brief yet illustrating introduction to the different concepts of algorithms, and [36j for a systematic 
introduction to Turing machines and their complexity. 
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When it comes to adding or multiplying rational numbers, there exist smart ways of designing 
the operations which can notoriously speed up the elementary algorithms, see for example [TB] . 
However, we will not search for the optimal upper bounds on the complexity of our algorithm, 
because our intention is just to prove that it is polynomial in certain quantities as claimed in 
Theorem |2] We just recall from [16] that k arithmetic operation^ (a.o. from now on) can be 
performed on rational inputs of bit length at most h, in time which is polynomial in k and h, that 
is in time {kh)^^^\ and the result of this sequence of a.o. is a rational number r e Q such that 
bl(r) < {kh)^^^\ 

Given a vector v G Q[i]'^, we define its bit length as 

bl(f ) = max{bl(aj), bl(6j) w = (ai + . . . , + \hk)}. 

7.1. Bit complexity of Computesign. Let h be an upper bound for the bit length of the input 
(a, 6, c, t, L, U) of Computesign. The algorithm performs a fixed number of arithmetic operations 
on the rational numbers which are its input. Hence, the bit complexity of Computesign is h^^^\ 

7.2. Bit complexity of LUquadratic. Let h be an upper bound for the bit length of the 
input {a,b,c, L,U) of LUquadratic. Until Step 11, LUquadratic performs a fixed number of 
arithmetic operations on the rational numbers which are its input (including two applications of 
Computesign). The bit complexity of LUquadratic until Step 11 is thus h'^^'^\ Each of the 
loops starting at Step 12 also performs a fixed number of arithmetic operations, but now the bit 
length of the number t2 invoked in COMPUTESIGN at line 19 grows with each loop. More precisely, 
after i iterations, 

bl(t2) < 0{t), 

and thus the maximum bit length in all the numbers appearing at the algorithm in the i-th loop 
is {h + The total bit complexity is thus 

^loops 

0{h) + 5Z (/i + if^''^ <ih + tlloops)^(^). 



From Remark 16 during an application of TrackSegment 

flloops < O logsmax 1, rii^ii n^^^ ii 

V V mm{||/t|| : < t < 1} 

Thus, the bit complexity of LUquadratic on inputs of bit length at most h is, during an appli- 
cation of TrackSegment, at most 

||/-(?||nrf3max{/i(/„Cf):0<t<l }^ ^ ^^'^ 
min{||/t|| : <t < 1} 



h + log2 max I 1, 



7.3. Bit complexity of ShortZero. Let h be an upper bound for the bit length of the input 
{z,e) of ShortZero. Steps 1 to 6 of ShortZero perform 0{n) a.o. on inputs of bit length 
bounded by h and thus these steps take time 

(n/i)0(i), 

which is also a bound for the bit length of x (in the notations of Algorithm ShortZero). The 
number of iterations the algorithm will perform is then at most 

0(log2||x||) <0(log2(n/i)). 



^By a.o. we mean an operation of the form +, — , x, /, or a comparison <, < or an assignment of a value to a 
variable, or computation of the integer part of a number. 
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at each step the bit length of a increases by a factor of 4, and checking the stopping criterion can 
be done in (n/i/log2(£))*^'-^-*. Hence the total bit complexity of the while loop is 

((loops 

^ {{nhl\og^{e)f^'^ ^0{i)) < (n/ilog2(£) + poops) (^(1) < {nh\o^^{e)f'^^\ 

i=\ 

Step 11 can then be done in (nh)'^^^\ Thus, the total bit complexity of ShortZero is (nh log2(£))*^*-^^ 

7.4. Bit complexity of TrackSegment. Let h be an upper bound for the bit length of the 
input (/, g, Zo) of TrackSegment. Let 5 > be the number of non-zero monomials in the dense 
representations of / and g. We assume that 

firnax = max{/i(/t, Ct) ^ < t < 1} < OO, 

which indeed implies that Co < oo and by Theorem [8] we know that TrackSegment actually 
produces an approximate zero of /. We now analyze the operations performed in each step of 
TrackSegment. 

(1) The operations before the while loop: 

• Steps 2,3,4: two squared-norm computations and one inner product computation. 
That is 0(5") a.o. with rationals of bit length max{/z, /} where / is an upper bound 
for the bit length of the multinomial coefficients (^') which appear in the definition 



of Bombieri-Weyl's product (see Section 2.1). Note that I < log{d\) < d^^^\ Thus 



max{/i, 1} < {h + (i)*^*^^^ < (hd)'^''^^ and the bit complexity of these steps is at most 

(Shd)'^^^\ The numbers they produce have bit length (Shd)'~'^^^ as well. 

Steps 1,5,6,7: a constant number of a.o. with rationals of bit length (Shd)'^^^^ is 

again (Shd)^^^^ (and the numbers produced have the bit length bounded by the same 

quantity). 



(2) Step 8 (number of loops): from Theorem |8] and Lemma 11, the number of loops is at most 
\79y/n + l(i^/^Co] , where Cq is the length of the path {ft,Ct) in the condition metric. For 
counting the bit complexity of each loop, let hi be h or the maximum bit length of the 
rational numbers Si,Zi,ti (whichever is greater), and let hmax = max{/ij} (we will prove 
latter that hmax < oo). The bit complexity of the i-th loop is bounded as follows. 

• Steps 9, 10, 11: a constant number of a.o. with rationals of bit length (hi)^^^^ is again 

(h^r^'i 

• Step 12: computation of the squared norm of a C"'^^ vector with rational coordinates 
of bit length bounded by h^: bit complexity (nhi)^^^^ and riy has bit length at most 
{nhif^^\ as well. 

• Step 13: computation of the derivative matrices of / and g at Zi, which is (nSdhi)'^^^^ 
using the elementary evaluation method (see [3] for a faster but more complicated 
one), and the bit length of the numbers is at most {nSdhi)^^^\ 

• Step 14: addition of two n x (n + 1) matrices with rational entries of bit length at 
most (nSdhi)'^^^^ is {nSdhi)'^^^\ then an inverse matrix computation is (nSdhi)'^^^^ 
using modular technique^ Indeed, computing of the inverse is equivalent to solving 
n + 1 systems of equations with rational coefficients. Each of these systems can be first 
normalized to systems with integer coefficients of size {nSdhi)'^^^\ which (according 
to, e.g., [18]) can be solved in time {nSdhi)^^^\ The total bit complexity of this step 
is thus {nSdhif^^\ 



6 

on the subject, as well as software and research articles. 



Exact linear algebra is a large research field, see 'http : //linalg . org/people . html for a list of people working 
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Step 15: 0(n^ log2((i)) arithmetic operations (the logg (i in this formula is needed to 
compute n^'"*"^ ^) with numbers of bit length {nSdhi)'^^^^ is again (nSdhi)'^^^\ and the 
bit length of a is again bounded by {nSdhi)^^^\ 

Steps 16, 17: computation of f{zi) and g{zi) is {nSdhi)'^^^^ because that is a bound 
for the bit length of the rational numbers appearing in the monomial expansion of /, g 
and also for the bit length of the coordinates of Zi. There are also a constant number 
of a.o. which is again {nSdhi)^^^\ 

Step 18: a constant number of a.o. is again {nSdhi)^^^\ 

Step 19: a matrix-vector product, O(n^) a.o. with rationals of bit length bounded by 
{nSdhi)0(^\ is again {nSdhifW^ 

Steps 20,21,22: a constant number of a.o. is again {nSdhi)'^^^\ 

Step 23: an application of LUquadr atic with input data whose bit length is bounded 
by {nSdhi)^^^\ according to Section 7.2 



n 



+ log2max ( 1, 



costs 

\\f - g\\nd^llrnax \ \ '^^^^ 



min{||/t|| :0<t< 1}, 
By hypothesis, the output of LUquadratic has the bit length bounded by /ij+i < 

• Step 24: a constant number of a.o. is again [nSdhmax]^^^'^ ■ 

• Step 25: a division of two rational numbers of the respective bit lengths (nSdhmax)'^^^^ 
and 

(7.1) bl(a) = bl(<^2) < 0{y/^\0g^fimax) 



costs (nSdhmax^Og^ fXmax)^''^^ ■ 

• Step 26: adding two n x [n + 1) matrices, O(n^) a.o. with rationals of bit length 

{nS dhmax ^Og2 fJ.max)'^^^^ haS bit complexity {nSdhmaxl^-maxT^^^ ■ 

• Step 27: as in Step 17, this takes time {nSdhi)'^^^\ 

• Step 28: solving a system of equations and adding two vectors with bit lengths bounded 
by {nS dhmax iog2 fimax)'^^^^ IS again {nSdhmaxfJ'max)^^^^ according to [18j. 

• Step 29: an application of ShortZe ro w ith input whose bit length is bounded by 
{nSdhmax log2 fJ'max)'^^^^ ■ From Section 7.3 , this has bit complexity {nSdhmax log2 l^max)'^^^'' ■ 

• Step 30, 31: a constant number of a.o. with rationals of bit length [nSdhmax log2 fJ'max)'^^^^ 

is {nSdhmax^Og2 Umax)^''^^ ■ 

(3) Step 33: One a.o. is again {nSdhmax^og2 fJ^max)'^^^^ ■ 
The bit complexity of TrackSegment is thus 

\\f - g\\nd^flrnax \\'^^^^ 



/ \\f - g\\nd^flmax \\ 

n'^ + logo 

max I 1, — : — , .. .. — , ) ) Co 

V mm{||/t|| : < t < Ijyy 



where h^ax is the maximum of h and the bit lengths of Si,ti and Zi. Now, all the Si and ti are 
numbers of the form m/2' where, from Remark 16 



bl(m) < bl(2') = l + l<0[ log2 max 1, 



min{||/4|| :0<t< 1}^ 

Thus, this is also an upper bound for the bit lengths of Sj and t^. As for that of Zi, note that from 
Lemma 15 we have that at each step i > 1, 



e 



h\{Zi) < O log2 < 0(log2(V^a)) < 0{\0g2{n^rnax))- 



1l3 
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Hence, we have 

h^ax <h+(^ (^log2 max (^1, o'^<^"< ) + ^og^iWmax^j ■ 

The bit complexity of TrackSegment is thus hnear in Cq and polynomial in the following quan- 
tities: 

• n, S, d, h, 

•log2(||/-^7||/min{||/,||:0<t<l}). 



8. Proof of Theorem [2] 



We first note that TrackSegment, performs the operations described by TrackSegment_Scheme, 
except for the use of Frobenius norm instead of operator norm in the computation of Xi,i- This 



follows directly from the description of the two algorithms and from lemmas 12 and 15 



Thus, from Theorem [8] and Lemma 11, TrackSegment has certified output. Moreover, its 
total bit complexity has been proved in section 7^ to satisfy the claim of Theorem [2j For the 
bound on the size of the output, let i = be the final step of the algorithm. Then, the output 
Zk+i of TrackSegment is the result of applying ShortZero to some {zk,e) where G Q[i 
and 

£^0 ^ Co ^ Ci 



n+1 



> 



> 



Co and Ci some constants. It follows from Lemma 15 that z^+i has integer coordinates of bit length 
at most 0(log2(n/i(/, ^fc+i))), as claimed. The proof is now complete. 



9. Experiments 

Our implementation of Algorithm |4] has been carried out in the top-level (interpreted) language 
of Macaulay2 [19j. The exact linear algebra routines and evaluation of polynomials are inherently 
slow and there are many engineering improvements that can be made to speed up the execution; 
yet the computation takes reasonable time on the examples of modest size. 

While more examples of computation along with the source code of the implementation are 
available at 

http : / / peo ple . math . gatech . edu/~aleykin3/RobustCHf7] 

here we describe two experiments. One of them involves a small family of equations, where most 
of the computation of the length of a homotopy path Co can be carried out by hand. The other 
comes from an application in enumerative geometry and showcases the class of problems that can 
benefit from the developed certified algorithms. 



9.1. Actual number of steps vs. condition length. In Lemma V\_ we claim that the number 
of steps (i.e. number of while loops) needed by Algorithm |4] is at most [79A/n -|- Id^/^Co]. In 
this section we consider a simple family of examples parametrized by m > where the value of 
Cq can be approximated by quadrature formulas and show how the bounds based on Cq compare 



to the actual performance of the algorithm. Note that from (2.3) the condition length of a path 



(with Ct given by a smooth curve of affine representatives) is 



Ml 
Wft 



^e{{ftjt)) , lie 



II/* 



+ 



IK* 



l(C*,C*)P 

iK*r 
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DftiCt] 

Ct* 



-1 



/Vd~i\\ft\\\\Ctr^-' 



\ 



/3;ii/tiiiictr"-' 



WQWJ 



In general, it is extremely hard to compute a priori Cq (even approximately). We consider here the 
simple case 

ft{xo, xi) = xl- {1 + mt)xl, (t = (1, Vl + mtY . 
Let s = 1 + mt. We can easily compute: 

- ' ' ft = -mxl; \\ftf = m^; {f^J^)=ms; 



\\ftr = l + s' 



Ct 



0, 



m 



IICt|P = i + ^; IK* 



m 
4^ 



l(4,0)P 



Thus, 



ii/*r 



llCtP 



\U\' 



1 



rriK 



+ 



1 



1 + S2 (1 + S2)2 4S(1 + S) 4(1 + S)2 

On the other hand. 



1 1 

+ 



(1+S2)2 4S(1+S)2' 



-1 



-2s 2v/^Y 
1 v^J 



-1 1 



2(l+s) 1+s 
1 



2v^(l+s) 1+s 



/ 2{l+s) 



2v^(l+s) 1+. 



V2s 



2s 



where to get the the last equality we compute the matrix norm by hand. With the change of 
variables s = 1 + mt we have then proved that 



Co(/t, Ci 



1 1 

+ 



1 + S2 2 4sl + S 



(is. 



It is not an easy task to find this integral exactly, but we can at least try to approximate with 
some quadrature formula. In Octave-produced Table [3] and Figure |3] we compare the values of 
upper and lower bounds 

LBound < tj(steps) < U Bound, where 

LBound = 28d^/^CoiftXt) ^ 79Co, 

U Bound = 79v/^^(i3/2Co(/t,Ct) =316Co(/i,Ct) 

for different choices of m > and the number of steps performed by our algorithm to follow the 
homotopy ft- 
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Table 3. Comparison of the bound of number of steps given by Lemma 11 and the 
actual number of steps in the example given hj ft = x1 — [1 + mt)xQ. 



m 




steps 


TTT3 


Uo/ steps 


iU 


Q 1 

oi 


1 S /I 

io4 


OO 1 


1 nc; 
i.yo 


on 


QS 
OO 


01 7 
/i / 


4o0 


o ni 


QD 
dU 


AO 


/O / 


/I sn 
4oU 


o nQ 
z.Uo 


An 


A K 


zOU 


Oiz 


o ni^ 
Z.UO 


ou 


4 ( 


ofin 


Oo / 


o nv 
z.U / 


fin 
du 


A n 
4y 


Ofin 


OOo 


o nfi 




c;n 


OVfi 


0/0 


o nfi 


sn 
oU 


oz 


OfiO 


c;nn 
OyU 


o nn 
Z.Uy 


nn 


Oo 


OSQ 

Zoo 


finQ 
dUo 


O 1 


1 nn 
iUU 


04 


ono 


fil e; 
oio 


oil 
Z.ii 


1 nnn 
iUUU 


11 
1 1 


Qn c; 
oyo 


S70 


o oi 
Z.Zi 


onnn 

ZUUU 


o4 


/I Ofi 
4Z0 


Q/1Q 

y4y 


0*? 

z.zo 


3000 


88 


446 


995 


2.23 


4000 


91 


457 


1027 


2.25 


5000 


93 


468 


1052 


2.25 


10000 


100 


499 


1129 


2.26 


20000 


106 


530 


1207 


2.28 


30000 


110 


547 


1252 


2.29 



9.2. An application to a problem in Schubert calculus. The computations of [28j confirmed 
the conjecture saying that the Galois group of a simple Schubert problem is the full symmetric 
group for "small" Grassmannians. These results produced using heuristic homotopy continuation 
methods take us far beyond the limitations of the symbolic methods. 

Table |4| a copy of [281 Table 1]), shows the number of solutions for the largest problem on 
G{k, n) and the number of permutations found in the Galois group by the algorithm sufficient to 
generate the full symmetric group. At the present all computations can be done within one day 
with a heuristic homotopy tracker employed. 



fc, n 


2,4 


2,5 


2,6 


2,7 


2,8 


2,9 


2,10 


solutions 


2 


5 


14 


42 


132 


429 


1430 


permutations 


4 


6 


5 


6 


7 


4 


7 



fc, n 


3,5 


3,6 


3,7 


3,8 


3,9 


4,6 


4,7 


4,8 


solutions 


5 


42 


462 


6006 


17589 


14 


462 


8580 


permutations 


4 


4 


5 


6 


7 


5 


5 


7 



Table 4. Galois group computation for simple Schubert problems in G{k,n). 



This problem falls naturally in the class where the certified algorithms of this paper can be 
applied. With the current implementation the algorithm of this paper can provide the status of a 
theorem to all of the computational results on up to Gr(2, 6): the cases that can be certified within 
a day appear in bold in Table |4j 

The corresponding runs of the algorithm for Gr(2, 6) involve tracking homotopies for six polyno- 
mial equations following the paths in and have input, output, and all intermediate approximate 
zeroes defined over Gaussian integers Z[i]. Due to the use of our Algorithm ShortZero to reduce 
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2.5 



U Bound /^{steps) 
LBound/^{steps) 




15000 20000 25000 30000 

m 



Figure 3. Comparison of the ratio between the actual number of steps and its 
lower and upper bound. 



the size of the integers in the intermediate steps, in this relatively large computation we do not 
encounter integers longer than six decimal digits amongst the coordinates of all approximate zeroes 
computed along all homotopy paths. 

Let us remark that the largest certifiable case is already beyond the reach of purely symbolic 
algorithms (the problem with 14 solutions in Gr(2,6) is characterized as "not computationally 
feasible" in |1Q]). There are several ways to push the frontier of provable results further. One is a 
low-level optimized implementation of our algorithm. Another is using a fast heuristic homotopy 
tracker to find the "interesting" paths (e.g., the ones that do not lead to a redundant permutation 
in the Galois group computation), break them up into a union of smaller pieces, and then execute 
a certified homotopy tracker for every small piece. The last step is trivially parallelizable and can 
be sped up in practice by using distributed computing. 
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10. Proof of Theorem [7] 

We recall first two lemmas [H Lemma 4 and Lemma 5]. The second of these two lemmas is 
recalled here in a less general version than the original. 

Lemma 17. Let ho, h E S, v E ^{d), zo,z G P(C"+"'^). Assume that Xi{ho,zo) < +oo. Assume 
moreover that 

a 

dR{zo, z) < 



ds{ho,h) < 



3d 



for some d < 1/V2. Then, 

^— < Xi{h, z) < ^ and 

1 + V2d ^ ^ ' ' 1 _ 

I, ^ (1 - 72 g)^ ^{ho,v,zo) 
Lp{ho,v,Zo)— — < Lp{h,v,z) < 



1 + V2d ~ ' ' ' ' ~ (1 - v^a)i+v^' 

Lemma 18. Let t^hsES,0<s<Tbea piece of a great circle in S, parametrized by arc-length. 
Let rjo G P(C"+"'^) be a projective zero of ho such that fi{ho,rio) < +oo. Assume that 

T < where <p = (p{ho, ho, r]o). 

Then, for < s < T , rjo can be continued to a zero rjs G P(C"^^) of hg in such a way that s ^ rjs 

is a C'l+^^p curve. Moreover, consider the curve s — {hg, hs,ris) , < s < T . Then, the following 
inequalities hold for every s G [0, T] : 

< <^{hg,hs,r]s) < 



ds{ho,hs) < log 



dVm ^l-d^/mx2iho,ho,Vo)s 

Now we proceed to the proof of Theorem jrj Recall that we have defined T = (^j|fj| > y/jj 
Consider the path 



(10.1) s K = jf^ cosis) + ^ sin(g), sG[0,T], 

' / g 
'11/11' llsli 

That is, hg is the arc-length parametrization of the short portion of the great circle joining g/\\g\ 
and //||/||. Note that ho = g as was defined in Theorem [?[ 
Let 

xi = xiigXo) = /i(^,Co),X2 = X2{g,gXo),0 = ^{9,9X0)- 

From (3.9) and Lemma 17 we get 

(1 - V2uo/2)^ 

(10.2) 0^ ' < V? < 



1 + V2uo/2 - (1 - V2«o/2)i+V2- 
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From (3.7) and (3.8), we have: 

(10.3) < fi{g, Co) 11^711 V^ + KaXoY < V2fi{g, Co)^ 

where for the last equahty we have used that ||/io|| = 1 and fi{g,C) ^ 1- The second inequahty of 
( |10.2D , together with ( |10.3D , then imphes ( |3.12D . 
It also follows that 

(10.4) 



T < 



< 



< 



pjot Prf^/V Pd3/^0 



Thus, Lemma 18 applies and we conclude that Co = ^70 can be continued to rjs E §(C"'+ ), a zero 
of hg, for < s < T. Now, note that hg is a reparametrization s = s{t) of the projection of 
/t = (1 — t)g + tf on §. That is, hg = ft{s)/\\ft{s)\\- Hence, Co can be continued to C* = Vs (t), a zero 
of ft as claimed in Theorem [7[ Note that Ci = Vs{i) = Vt- Moreover, Lemma 18 and (3.7) also 
imply that, for < s < T, 



(10.5) 
and that 



< K^s, Vs) II {hs, r]s) ||t(;,^,,^)§xP(C"+i) 



< 



1 -P(i3/2^S' 



dniCoXi) = dR{r]o,r]T) < 



1 



(10.6) 



V2d^/'Xi 



V2dy^xi 

1 - (1 - c')^/^ 



Pd^/-^ip J 



< 

i |Ta2l i,p?2li 



We have seen (10.4) that T < c'{Pd^^^(f>) ^. Now, (p = X1X2 and X2 > ||^s|| = 1? which implies 

3a 



< 



Pd^/^Xi Pdy^Kg, Co) S 2V2dy^fL{g, Co 



T < 

Note that this last inequality, (10.6) and Lemma [17| imply 

Kg, Co) 



(10.7) 

Thus, 
(10.8) 
and hence 



Kg, Co) 



</i(/,Ci)< 



1 — a 



1 + a 

KgXo) > (1 - a)KfXi 



(10.9) 
Note that 



< 



Pc?3/2(l-a)^(/,Ci) PI 2f/3/2^(/,Cl)' 



T < 



0^^(2:0, Ci)Kf, Ci) < (c^ij(2o, Co) + dniCo, Ci))Kf, Ci] 

Uq a 



< 



+ 



3.9M10.6) \2d^/^KgXo) V2dy^fi{g,Co) 

1 



M/,Ci) 



< 

JTotI 



+ 
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Our choice of a is such that the right-hand term in this last equation is at most 5uo/(f^'^. Hence, 
we have 

Suo 



(10.10) 



From (10.10) and (10.9), Lemma 17 then yields 



(10.11) 



1 + a/25mo 



<Xi< 



1 - y/26uo ■ 



Moreover, from Lemma [6| (10.10) implies that zq is an approximate zero of / with associated zero 
(i. In particular, 

6uo 



(10.12) 



ci^(iVp(/)(zo),Ci)<^^%^< 



From this and (3.14) we have 

dR{S,(i) < dRCz,N^{f){zo)) + dR{N^{f){zo),(i) < 
(1 - 5)uo , 5uo 



(1 - 6)uo , Suq 



+ 



< 



< 



2o?3/>(/,Ci) 



((1-5) + 5) 



2ci3/V(/,Ci)' 



proving (3.15). 



As for (3.13), first note that from Lemma 19 below. 



(10.13) 



H{hs, r]s) II {hs, Vs) ||T(h^,,^)§xP(C"+i) ds. 



Now, this last equality and (10.5) imply: 

Co(/t,Ct)> 



1 + Pdy^<fs 



ds 



ln(l + Pd^/^^T) _ ln(l + Pd^/'^^T) 
Because ln(l + t)/t is a decreasing function of t > and 



( [Ta2l i 



we conclude that 
(10.14) 



;i-y2Mo/2)^ ^ {l-V2uo/2)^ 



CoiftXt) > > 



<c' 



;i- V2Mo/2)i+^ln(l + c') 



On the other hand, using again (10.13), we have 
(10.15) CoiftXt) < 



1-Prf3/2(^S 



ds 
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Note that (10.14) and (10.15) prove ( |3.13[ ). This finishes the proof of Theorem [7| 



We have to prove a lemma that has been used in the proof of Theorem [7} and which is nothing 
but a change of variables: 

Lemma 19. In the notation of the proof of Theorem^ we have: 

Co{ft,(t)= / /i(/ls,?/s)||(/i.,?7s)||T(^^,^)§xP(C"+i)C?S- 

Jo 

Proof. One can just apply the change of variables formula to the change of variables s = s{t) (so 
that hs(^t) = /t/ II /ill ris(t) = Ct) and, after a long computation prove that the two integrals 
of the lemma are equal. However, we prefer the following geometric argument. The quantity 
Co(/f,Ci) is by definition the length of the path (/t/||/t||, Ct) when § x P(C"'+^) is endowed with 
the condition metric, resulting from multiplying the usual product metric by the square of the 
condition number at each pair {f,z). Now, as a length, it is independent of the parametrization 
and thus Co{ft, Ct) = Co{hs,ris). This is exactly the claim of the lemma. □ 
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